DEMO: H wavefunctions

Requirements | Importing libraries

import matplotlib, matplotlib.pyplot as plt
import seaborn as sns
import numpy as np

from scipy.constants import physical_constants
from matplotlib import cm, colors
import plotly.graph_objects as go
from mpl_toolkits.mplot3d import Axes3D

%matplotlib inline
%config InlineBackend.figure_format = 'retina'

try:
    from google.colab import output
    output.enable_custom_widget_manager()
    print('All good to go')
except:
    print('Okay we are not in Colab just proceed as if nothing happened')

Okay we are not in Colab just proceed as if nothing happened

1. Describing a radial function \(R_{nl}(r)\)

\[R_{nl} = \sqrt{ \Big( \frac{2}{n a_0} \Big)^3 \frac{(n-l-1)!}{2n(n+l)!}} \cdot e^{-\frac{r}{n a_0}} \cdot \Big( \frac{2 r}{n a_0} \Big)^l \cdot L^{2l+1}_{n-l-1}\Big( \frac{2 r}{n a_0} \Big)\]

def radial_function(r, 
                    n=1, 
                    l=0):
    '''Rnl(r) normalized radial function
    r: radial distance Float (0, inf)
    n: principal quantum number Int (1,2,3... inf)
    l: angular quantum number Int (0,1,2,... n-1)
    '''
    from scipy.special import  genlaguerre
    
    
    a0 = 1 # Bohr radius equal to 5.29e-11 m

    prefactor = np.sqrt( ((2 / n * a0) ** 3 * (np.math.factorial(n - l - 1))) / (2 * n * (np.math.factorial(n + l))) )

    laguerre = genlaguerre(n - l - 1, 2 * l + 1)

    p = 2 * r / (n * a0)
    

    return  prefactor * np.exp(-p / 2) * (p ** l) * laguerre(p)

r = np.linspace(0, 20, 1000)

Rnl = radial_function(r, n=2, l=0)

#plt.plot(r, Rnl)
#plt.plot(r, Rnl**2)
plt.plot(r, r**2 * Rnl**2)
plt.xlabel(r'$r [a_0]$')
plt.ylabel(r'$R_nl(r) r^2$')

/tmp/ipykernel_15473/1039204017.py:14: DeprecationWarning: `np.math` is a deprecated alias for the standard library `math` module (Deprecated Numpy 1.25). Replace usages of `np.math` with `math`
  prefactor = np.sqrt( ((2 / n * a0) ** 3 * (np.math.factorial(n - l - 1))) / (2 * n * (np.math.factorial(n + l))) )

Text(0, 0.5, '$R_nl(r) r^2$')

2. Describing an angular function | Spherical harmonic \(Y_{l}^{m}(\theta, \varphi)\)

• We will make use of sph_harm function from Scipy

• We can also build up spherical harmonics using Associated Legendre function using lpvm

\[ Y_{lm}(\phi,\theta) = \sqrt{\frac{2l+1}{4\pi} \frac{(l-m)!}{(l+m)!} } P_{lm}(cos \phi) \cdot e^{im\theta} \]

from scipy.special import sph_harm

sph_harm(0, 0, np.pi, np.pi) # test a few spherical harmonics

(0.28209479177387814+0j)

def plot_Yml(l, m):
    '''Visualizing spherical harmonics using sph_harm funcion from Scipy special function library
     Note that the name of angles is different from the notation adopted in QM textbooks!
     l: angular quantum number (0,1,2,...)
     m: magnetic quantum number (-l,...l)
     '''

    #Creade grid of phi and theta angles for ploting surface mesh
    phi, theta = np.linspace(0, np.pi, 100), np.linspace(0, 2*np.pi, 100)
    phi, theta = np.meshgrid(phi, theta)

    #Calcualte spherical harmonic with given m and l
    Ylm = sph_harm(m, l, theta, phi) 
    R=abs(Ylm)
    
    # Let's normalize color scale
    fcolors    = Ylm.real
    fmax, fmin = fcolors.max(), fcolors.min()
    fcolors    = (fcolors - fmin)/(fmax - fmin)

    # Since we want to plot on cartesian reference frame we will use cartesian coordiniates x, y, z using R as the absolute value of Yml
    # Try a plot without R part. 
    fig = go.Figure(data=[go.Surface(x=R*np.sin(phi) * np.cos(theta), 
                                     y=R*np.sin(phi) * np.sin(theta), 
                                     z=R*np.cos(phi), 
                                     surfacecolor=fcolors,
                                     colorscale='balance')])

    # Show the plot
    fig.update_layout(title=fr'$Y_{l, m}$', autosize=False,
                      width=700, height=700,
                      margin=dict(l=65, r=50, b=65, t=90)
    )
    fig.show()

plot_Yml(l=1, m=0)

3. Describing the normalized probability as wavefunction squared \(|\psi _{nlm}(r,\theta ,\varphi)|^2\)

\[\psi_{nlm} = R_{nl}(r) \cdot Y_{l}^{m}(\theta, \varphi)\]

def normalized_wavefunction(n, 
                            l, 
                            m, 
                            max_r = 10):

    '''Ψnlm(r,θ,φ) normalized wavefunction
    by definition of quantum numbers n, l, m and a bohr radius augmentation coefficient'''

    # Set coordinates grid to assign a certain probability to each point (x, y) in the plane
    x = y = np.linspace(-max_r, max_r, 501)
    x, y = np.meshgrid(x, y)
    r = np.sqrt((x ** 2 + y ** 2))

    # Ψnlm(r,θ,φ) = Rnl(r).Ylm(θ,φ)
    psi = radial_function(r, n, l) * sph_harm(m, l, 0, np.arctan(x / (y + 1e-10)))

    return np.abs(psi) ** 2

4. Plotting wavefunction electron probability density plots

Light shaded areas in the orbital cross-sections represent a high probability of a particle being present in that region.

def plot_orbital_2d(n,l,m, max_r = 10):
    fig, ax = plt.subplots()
  
    im = ax.contour(normalized_wavefunction(n, l, m, max_r=max_r), 16, cmap='rocket', extent=[-max_r, max_r,-max_r, max_r])
    ax.set_title(r'$|\psi_{{({0}, {1}, {2})}}|^2$'.format(n, l, m), fontsize=15)
    ax.set_aspect('equal')

plot_orbital_2d(5,3,0, max_r=20)

/tmp/ipykernel_15473/1039204017.py:14: DeprecationWarning:

`np.math` is a deprecated alias for the standard library `math` module (Deprecated Numpy 1.25). Replace usages of `np.math` with `math`